%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     This program simulates      %
%       MSE for Figure 2          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Set up 
addpath 'A:\VAM2\cpatel\code\ado';              %Set to folder given by Stata global "ado"
clear all; 
rng(40)

city='NYCms';                                   %Select city ("NYC" or "NYCms")
vtype='homo';

ols=5; 	  %Number of VAMs
S=500;    %Number of simulations
fill = 1; %Load inputs and save mat file

%%%%%%%%%% LOAD DATA %%%%%%%%%%%
if strcmp(city, 'NYCms')
    dir=['D:\nyc\VAM2\cpatel\build\'];          %Set to folder given by Stata global "code_build_NYCms"
    outdir='A:\VAM2\cpatel\output\tables\';     %Set to folder given by Stata global "tables"
    fg = '6';
    lg = '6';
    fy = '2016';
    ly = '2018';
    bw='_ccftcovs_bwall';
    ptype='form2';
end
if strcmp(city, 'NYC')
    dir=['D:\nyc\VAM2\cpatel\build\'];          %Set to folder given by Stata global "code_build_NYC"
    outdir='A:\VAM2\cpatel\output\tables\';     %Set to folder given by Stata global "tables"
    fg = '9';
    lg = '9';
    fy = '2012';
    ly = '2014';
    bw='_ccftcovs_bwall';
    ptype='form2';
end
if fill == 1
    samp = [fg '_' lg '_' fy '_' ly '_e1_q1'];
    sects=dlmread([dir city '_sectors_' ptype '_' samp bw '.csv']);
    J = size(sects,1);
    J1 = J+1;
    L = sum(sects(:,1));
    dim = J+L+J*L;

    data_temp = [];
    V_raw = sparse(dim,dim);
    for i=1:J+2          
         diagdiag = 'diagdiag_m';
         b_block = dlmread([dir diagdiag '\' city '_' ptype bw '_b' num2str(i) '.csv']);
         V_block = dlmread([dir diagdiag '\' city '_' ptype bw '_V' num2str(i) vtype '.csv']);
         V_block_diag = diag(V_block);
         V_block = diag(V_block_diag);
         data_temp = [data_temp , b_block];
         if i==1
            V_raw(1:J,1:J) = V_block;
         end
         if i>1
            V_raw(J+L*(i-2)+1:J+L*(i-1),J+L*(i-2)+1:J+L*(i-1)) = V_block;
         end
         i
     end
    
        %OLS VAM alphas;
        alpha1 = dlmread([dir city  '_b_unc.csv'])';
        alpha2 = dlmread([dir city  '_b_conv.csv'])';
        alpha3 = dlmread([dir city  '_b_risk.csv'])';
        alpha4 = dlmread([dir city  '_b_rcvam.csv'])';
        alpha5 = dlmread([dir city  '_b_conv_3yrlag.csv'])';

        %OLS VAM variances;
        V_alpha11 = dlmread([dir city  '_V_unc.csv'])';
        V_alpha22 = dlmread([dir city  '_V_conv.csv'])';
        V_alpha33 = dlmread([dir city  '_V_risk.csv'])';
        V_alpha44 = dlmread([dir city  '_V_rcvam.csv'])';
        V_alpha55 = dlmread([dir city  '_V_conv_3yrlag.csv'])';
        
        %OLS VAM covariances;
        V_alpha14   = dlmread([dir city '_' ptype '_alpha_unc_rcvam_sampling_cov_m.csv']);
        V_alpha24   = dlmread([dir city '_' ptype '_alpha_conv_rcvam_sampling_cov_m.csv']);
        V_alpha34   = dlmread([dir city '_' ptype '_alpha_risk_rcvam_sampling_cov_m.csv']);
        V_alpha12   = dlmread([dir city '_' ptype '_alpha_unc_conv_sampling_cov_m.csv']);
        V_alpha32   = dlmread([dir city '_' ptype '_alpha_risk_conv_sampling_cov_m.csv']);
        V_alpha13   = dlmread([dir city '_' ptype '_alpha_unc_risk_sampling_cov_m.csv']);   
        V_alpha15   = dlmread([dir city '_' ptype '_alpha_unc_conv_3yrlag_sampling_cov_m.csv']);
        V_alpha25   = dlmread([dir city '_' ptype '_alpha_conv_conv_3yrlag_sampling_cov_m.csv']);
        V_alpha35   = dlmread([dir city '_' ptype '_alpha_risk_conv_3yrlag_sampling_cov_m.csv']);
        V_alpha45   = dlmread([dir city '_' ptype '_alpha_rcvam_conv_3yrlag_sampling_cov_m.csv']);
        V_alpha16   = dlmread([dir city '_' ptype '_alpha_unc_conv_3yrlag_sampling_cov_m.csv']);

        %RF variance and RF-OLS covariances;
        V_rho        = dlmread([dir city '_rf_V_1.csv']); 
        V_rho_alpha1 = dlmread([dir city '_V_alpharho_unc_JL_new_m.csv'])';  
        V_rho_alpha2 = dlmread([dir city '_V_alpharho_conv_JL_new_m.csv'])';
        V_rho_alpha3 = dlmread([dir city '_V_alpharho_risk_JL_new_m.csv'])';
        V_rho_alpha4 = dlmread([dir city '_V_alpharho_rcvam_JL_new_m.csv'])';
        V_rho_alpha5 = dlmread([dir city '_V_alpharho_conv_3yrlag_JL_new_m.csv'])';
    
    %Lottery weights;
    ZZ=csvread([dir city '_' ptype bw '_ZZ.csv']);
    inv_ZZ=inv(ZZ);

    %Hyperparameters;
    hyperp=csvread([dir city '_hyperparameters_7.csv']);
    mu=hyperp(1:6,1);
    sigma=hyperp(7,1);
    mu=mu([6 1 2 3 5 4],:);
    
    save([dir city '_' ptype '_matlab_filled.mat'])
else
    load([city '_' ptype '_matlab_filled.mat'])
end

%Format;
V_short=V_raw(1:J+L+J*L,1:J+L+J*L);
rho_hat=data_temp((J+1):(J+L))';
Pi_raw=data_temp(J+L+1:J+L+L*J);
Pi_raw=reshape(Pi_raw,[L,J]);
Pi=[Pi_raw,-sum(Pi_raw,2)];
V = V_short(J+1:J+L,J+1:J+L);

%Parse sectors;
if strcmp(city, 'NYC')
    screened = sects(:,2);
    sectors = [screened;0];
end
if strcmp(city, 'NYCms')
    screened = sects(:,2);
    sectors = [screened;0];
end
T=size(sectors,2);

%Remove schools not represented in older lag sample;
if strcmp(city,'NYC') 
    A = [477 468 451 376 331 325 322 183 156 71 5];
    for i = 1:length(A)
     Pi(:,A(i))=[];
     sectors(A(i),:)=[];
    end
    A = [376 331 325 322 183 156 71 5];
    for i = 1:length(A)
     Pi(A(i),:)=[];
     V_rho(A(i),:) = [];
     V_rho(:,A(i)) = [];
    end
    J=J-11;
    J1=J1-11;
    L=L-8;
end    
if strcmp(city,'NYCms')
    A = [590 522];
    for i = 1:length(A)
     Pi(:,A(i))=[];
     sectors(A(i),:)=[];
    end
    J=J-2;
    J1=J1-2;
end  

%Joint OLS covariance matrix;
V_alpha_full = [V_alpha11  V_alpha12  V_alpha13  V_alpha14  V_alpha15; ...
                V_alpha12' V_alpha22  V_alpha32' V_alpha24  V_alpha25; ...
                V_alpha13' V_alpha32  V_alpha33  V_alpha34  V_alpha35; ...
                V_alpha14' V_alpha24' V_alpha34' V_alpha44  V_alpha45; ...
                V_alpha15' V_alpha25' V_alpha35' V_alpha45' V_alpha55]; 
V_alpha_full = nearestSPD(V_alpha_full);

%Joint RF-OLS covariance matrix;
Sigma = [V_alpha11  V_alpha12  V_alpha13  V_alpha14  V_alpha15  V_rho_alpha1'; ...
         V_alpha12' V_alpha22  V_alpha32' V_alpha24  V_alpha25  V_rho_alpha2'; ...
         V_alpha13' V_alpha32  V_alpha33  V_alpha34  V_alpha35  V_rho_alpha3'; ...
         V_alpha14' V_alpha24' V_alpha34' V_alpha44  V_alpha45  V_rho_alpha4'; ...
         V_alpha15' V_alpha25' V_alpha35' V_alpha45' V_alpha55  V_rho_alpha5'; ...
         V_rho_alpha1 V_rho_alpha2 V_rho_alpha3 V_rho_alpha4 V_rho_alpha5 V_rho]; 
Sigma_rho_alpha = {V_rho_alpha1 V_rho_alpha2 V_rho_alpha3 V_rho_alpha4 V_rho_alpha5};
Sigma = nearestSPD(Sigma);   

%Parse OLS-VAMs;
alphas = {alpha1,alpha2,alpha3,alpha4,alpha5};
V_alphas = {V_alpha11,V_alpha22,V_alpha33,V_alpha44,V_alpha55};     
X=[sectors,alpha1,alpha2,alpha3,alpha4,alpha5];
K=5;

%%%%% Sigma alpha %%%%%%
w=1;
C = zeros(5,5);
C(1,1) = var(alpha1,w) - mean(diag(V_alpha11));
C(2,2) = var(alpha2,w) - mean(diag(V_alpha22));
C(3,3) = var(alpha3,w) - mean(diag(V_alpha33));
C(4,4) = var(alpha4,w) - mean(diag(V_alpha44));    
C(5,5) = var(alpha5,w) - mean(diag(V_alpha55));    
cov12 = cov(alpha1,alpha2);
cov13 = cov(alpha1,alpha3);
cov14 = cov(alpha1,alpha4);
cov15 = cov(alpha1,alpha5);
cov23 = cov(alpha2,alpha3);
cov24 = cov(alpha2,alpha4);
cov25 = cov(alpha2,alpha5);
cov34 = cov(alpha3,alpha4);
cov35 = cov(alpha3,alpha5);
cov45 = cov(alpha4,alpha5);
C(1,2) = cov12(1,2) - mean(diag(V_alpha12));
C(1,3) = cov13(1,2) - mean(diag(V_alpha13));
C(1,4) = cov14(1,2) - mean(diag(V_alpha14));
C(1,5) = cov15(1,2) - mean(diag(V_alpha15));
C(2,3) = cov23(1,2) - mean(diag(V_alpha32));
C(2,4) = cov24(1,2) - mean(diag(V_alpha24));
C(2,5) = cov25(1,2) - mean(diag(V_alpha25));
C(3,4) = cov34(1,2) - mean(diag(V_alpha34));
C(3,5) = cov35(1,2) - mean(diag(V_alpha35));
C(4,5) = cov45(1,2) - mean(diag(V_alpha45));
for c=1:5
    for d=1:5
        if c ~= d && c > d
            C(c,d) = C(d,c);
        end
    end
end
C
Sigma_alpha = C;
Sigma_alpha = nearestSPD(Sigma_alpha)

%%%%%% SIMULATE %%%%%%%

%Initialize;
sim_Pi=zeros(L,J+1,S);
sim_beta=zeros(J+1,S);
sim_rho=zeros(L,S);
sim_mu=zeros(S,size(mu,1));
sim_lambda=zeros(S,1);
sim_sigma=zeros(S,1);
sim_sigma_beta=zeros(S,1);
sim_post_alpha=zeros(J+1,K,S);
sim_post_hybrid1=zeros(J+1,S);
sim_post_hybrid2=zeros(J+1,S);
sim_post_hybrid3=zeros(J+1,S);
sim_post_hybrid4=zeros(J+1,S);
sim_post_hybrid5=zeros(J+1,S);
sim_post_hybrid6=zeros(J+1,S);
sim_post_hybrid7=zeros(J+1,S);
sim_post_hybrid8=zeros(J+1,S);

%Draw alphas;
sim_alpha = zeros(J+1,K,S);
for s=1:S
    sim_alpha(:,:,s)=mvnrnd(zeros(1,K),Sigma_alpha,J+1); 
end

%Store OLS;
sim_X_star=repmat(X,1,1,S);
sim_X_star(:,T+1:T+K,:)=sim_alpha;

%Redraw errors to get OLS posterior variance;
for sim=1:100
    sim_eps = mvnrnd(zeros(1,K*J1+L),Sigma,S)';
    sim_eps = sim_eps(1:ols*J1,:);
    sim_eps = reshape(sim_eps,[J+1,K,S]);
    for s=1:S
        sim_alphahat_var(:,:,s) = sim_alpha(:,:,sim) +sim_eps(:,:,s);
    end
        
    %%%demean
    sim_alphahat_var=sim_alphahat_var-repmat(mean(sim_alphahat_var),J+1,1);
    
    %%% shrink alpha_hats
    for k=1:K
       lambda(:,k)=Sigma_alpha(k,k)./(Sigma_alpha(k,k)+diag(V_alphas{k}));
       sim_post_alpha_var(:,k,:)=lambda(:,k).*sim_alphahat_var(:,k,:);
    end
    ols_var(sim,:)=mean(var(sim_post_alpha_var,1,3));
end
ols_var=mean(ols_var);

%Redraw errors to get IV posterior variances; 
for sim=1:100
    sim_eps = mvnrnd(zeros(1,K*J1+L),Sigma,S)';
    sim_e = sim_eps(ols*J1+1:ols*J1+L,:);
    sim_eps = sim_eps(1:ols*J1,:);
    sim_eps = reshape(sim_eps,[J+1,K,S]);
    sim_nu=sigma*randn(J+1,S);
     
    for s=1:S
        sim_beta(:,s)=sim_X_star(:,:,sim)*mu+sim_nu(:,sim);
        sim_rho_var(:,s)=Pi*sim_beta(:,s)+sim_e(:,s);  
        sim_alphahat_var(:,:,s) = sim_alpha(:,:,sim) +sim_eps(:,:,s);
        for v=1:ols
                 %estimate hyperparameters
                 sim_mu(s)=((sim_alphahat_var(:,v,s)'*Pi'*inv_ZZ*Pi*sim_alphahat_var(:,v,s))\(sim_alphahat_var(:,v,s)'*Pi'*inv_ZZ*sim_rho_var(:,s)))';
                 res=sim_rho_var(:,s)-Pi*sim_alphahat_var(:,v,s)*sim_mu(s)';
                 sim_sigma(s)=sqrt(res'*ZZ*res/trace(Pi'*ZZ*Pi));
                 sim_sigma_beta(s)=sqrt(var(sim_beta(:,s)));
                 % hybrid posteriors (estimated hyperparameters)
                 sim_sigma_alpha = sqrt(var(sim_alphahat_var(:,v,s)) - mean(diag(V_alphas{v})));
                 V_beta_beta=((sim_mu(s)^2)*(sim_sigma_alpha^2)+(sim_sigma(s)^2))*eye(J1);
                 V_beta_rho=((sim_mu(s)^2)*(sim_sigma_alpha^2)+(sim_sigma(s)^2))*(Pi');
                 V_beta_alpha=sim_mu(s)*(sim_sigma_alpha^2)*eye(J1);   
                 V_rho_rho=((sim_mu(s)^2)*(sim_sigma_alpha^2)+(sim_sigma(s)^2))*(Pi*(Pi'))+V_rho;
                 V_rho_alpha=(sim_mu(s)*(sim_sigma_alpha^2))*Pi + Sigma_rho_alpha{v};
                 V_alpha_alpha=(sim_sigma_alpha^2)*eye(J1) + V_alphas{v};
                 Omega_rho = (V_beta_rho - V_beta_alpha*inv(V_alpha_alpha)*(V_rho_alpha'))*inv(V_rho_rho - V_rho_alpha*inv(V_alpha_alpha)*(V_rho_alpha'));
                 Omega_alpha=(V_beta_alpha - V_beta_rho*inv(V_rho_rho)*V_rho_alpha)*inv(V_alpha_alpha - (V_rho_alpha')*inv(V_rho_rho)*V_rho_alpha);

                 sim_post_hybrid_var(:,v,s)=Omega_rho*sim_rho_var(:,s) + Omega_alpha*sim_alphahat_var(:,v,s);   
        end
    s    
    end
    iv_var(sim,:)=mean(var(sim_post_hybrid_var,1,3));
end
iv_var=mean(iv_var);

%Jointly draw OLS and RF errors;
sim_eps = mvnrnd(zeros(1,K*J1+L),Sigma,S)';
sim_e = sim_eps(ols*J1+1:ols*J1+L,:);
sim_eps = sim_eps(1:ols*J1,:);
sim_eps = reshape(sim_eps,[J+1,K,S]);

%Form alpha hats;
sim_alphahat=sim_alpha+sim_eps;
sim_alphahat=sim_alphahat-repmat(mean(sim_alphahat),J+1,1);

%Store OLS;
sim_X=repmat(X,1,1,S);
sim_X(:,T+1:T+K,:)=sim_alphahat;

%Draw forecast residual;
sim_nu=sigma*randn(J+1,S);

%Form beta and RF (hold Pi fixed);
for s=1:S
   sim_beta(:,s)=sim_X_star(:,:,s)*mu+sim_nu(:,s);
   sim_Pi(:,:,s)=Pi; 
   sim_rho(:,s)=sim_Pi(:,:,s)*sim_beta(:,s)+sim_e(:,s); 
end

%Shrink alpha_hats;
for k=1:K
   lambda(:,k)=Sigma_alpha(k,k)./(Sigma_alpha(k,k)+diag(V_alphas{k}));
   sim_post_alpha(:,k,:)=lambda(:,k).*sim_alphahat(:,k,:);
end

%%% calculate posteriors using VAMs as mediators one at a time 
for v=2:6 
    sim_mu=zeros(S,1);
    for s=1:S    
        %estimate hyperparameters
        sim_mu(s)=((sim_X(:,v,s)'*sim_Pi(:,:,s)'*inv_ZZ*sim_Pi(:,:,s)*sim_X(:,v,s))\(sim_X(:,v,s)'*sim_Pi(:,:,s)'*inv_ZZ*sim_rho(:,s)))';
        res=sim_rho(:,s)-sim_Pi(:,:,s)*sim_X(:,v,s)*sim_mu(s)';
        sim_sigma(s)=sqrt(res'*ZZ*res/trace(sim_Pi(:,:,s)'*ZZ*sim_Pi(:,:,s)));
        sim_sigma_beta(s)=sqrt(var(sim_beta(:,s)));

        % hybrid posteriors (estimated hyperparameters)
        sim_sigma_alpha = sqrt(var(sim_X(:,v,s)) - mean(diag(V_alphas{v-1})));
        V_beta_beta=((sim_mu(s)^2)*(sim_sigma_alpha^2)+(sim_sigma(s)^2))*eye(J1);
        V_beta_rho=((sim_mu(s)^2)*(sim_sigma_alpha^2)+(sim_sigma(s)^2))*(sim_Pi(:,:,s)');
        V_beta_alpha=sim_mu(s)*(sim_sigma_alpha^2)*eye(J1);   
        V_rho_rho=((sim_mu(s)^2)*(sim_sigma_alpha^2)+(sim_sigma(s)^2))*(sim_Pi(:,:,s)*(sim_Pi(:,:,s)'))+V_rho;
        V_rho_alpha=(sim_mu(s)*(sim_sigma_alpha^2))*sim_Pi(:,:,s) + Sigma_rho_alpha{v-1};
        V_alpha_alpha=(sim_sigma_alpha^2)*eye(J1) + V_alphas{v-1};
        Omega_rho = (V_beta_rho - V_beta_alpha*inv(V_alpha_alpha)*(V_rho_alpha'))*inv(V_rho_rho - V_rho_alpha*inv(V_alpha_alpha)*(V_rho_alpha'));
        Omega_alpha=(V_beta_alpha - V_beta_rho*inv(V_rho_rho)*V_rho_alpha)*inv(V_alpha_alpha - (V_rho_alpha')*inv(V_rho_rho)*V_rho_alpha);

        if v==2
           sim_post_hybrid3(:,s)=Omega_rho*sim_rho(:,s) + Omega_alpha*sim_X(:,v,s);
        end
        if v==3 
           sim_post_hybrid4(:,s)=Omega_rho*sim_rho(:,s) + Omega_alpha*sim_X(:,v,s);
        end
        if v==4
           sim_post_hybrid5(:,s)=Omega_rho*sim_rho(:,s) + Omega_alpha*sim_X(:,v,s);
        end
        if v==5 
           sim_post_hybrid6(:,s)=Omega_rho*sim_rho(:,s) + Omega_alpha*sim_X(:,v,s);
        end
        if v==6
           sim_post_hybrid7(:,s)=Omega_rho*sim_rho(:,s) + Omega_alpha*sim_X(:,v,s);
        end
    end
end

% put quality measures into a cell
tab=zeros(3,4,S);
measure{1} = sim_beta;
for k=1:K
   measure{k+1} = sim_post_alpha(:,k,:);
end
measure{K+2}=sim_post_hybrid3; 
measure{K+3}=sim_post_hybrid4; 
measure{K+4}=sim_post_hybrid5; 
measure{K+5}=sim_post_hybrid6;
measure{K+6}=sim_post_hybrid7;
measure{K+7}=sim_post_hybrid8;

%%% simulate closures
for s=1:S
   beta=measure{1}(:,s);
   
   for k=1:size(measure,2)
       j=1;

       %correlation with truth
       tab(k,j,s)=corr(measure{k}(:,s),beta);
       tab(k,j+1,s)=corr(measure{k}(:,s),beta,'Type','Spearman');
       reliability=fitlm(measure{k}(:,s),beta);
       tab(k,j+2,s)= reliability.Coefficients.Estimate(2);
       j=j+3;

       %closing lowest ranked school
       [~,l]=min(measure{k}(:,s));
       tab(k,j,s)=mean(beta)-beta(l);
       m=median(measure{k}(:,s));
       d=measure{k}(:,s)>m;
       tab(k,j+1,s)=mean(beta.*d)/mean(d)-beta(l);
       m=prctile(measure{k}(:,s),80);
       d=measure{k}(:,s)>m;
       tab(k,j+2,s)=mean(beta.*d)/mean(d)-beta(l);
       j=j+4;

       %closing lowest in district
       j=j+4;

       %closing bottom 5%
       m=prctile(measure{k}(:,s),5);
       l=measure{k}(:,s)<=m;
       tab(k,j,s)=mean(beta)-mean(beta(l));
       m=median(measure{k}(:,s));
       d=measure{k}(:,s)>m;
       tab(k,j+1,s)=mean(beta.*d)/mean(d)-mean(beta(l));
       m=prctile(measure{k}(:,s),80);
       d=measure{k}(:,s)>m;
       tab(k,j+2,s)=mean(beta.*d)/mean(d)-mean(beta(l));
       j=j+4;

       %closing bottom 5% in district
       j=j+4;
             
        %Close bottom 10;
        m=prctile(measure{k}(:,s),10);
        l=measure{k}(:,s)<=m;
        tab(k,j,s)=mean(beta)-mean(beta(l));
        j=j+1;
        
        %residual;
        reliability=fitlm(measure{k}(:,s),beta);
        tab(k,j,s)= std(reliability.Residuals.Raw);
        j=j+1;

        %MSE total, bias, variance;
 		tab(k,j,s) = mean((measure{k}(:,s)-beta).^2);
   end
end

% store variance components from first loop
variance = [ols_var(1,1) ols_var(1,1) ols_var(1,2) ols_var(1,3) ...
           ols_var(1,4) ols_var(1,5) ols_var(1,5) ols_var(1,1) ...
           ols_var(1,1) ols_var(1,1) ols_var(1,1) ols_var(1,1); ...
           iv_var(1,1) iv_var(1,1) iv_var(1,2) iv_var(1,3) ...
           iv_var(1,4) iv_var(1,5) iv_var(1,1) iv_var(1,2) ...
           iv_var(1,3) iv_var(1,4) iv_var(1,5) iv_var(1,5)];

% store
for s=1:S
   beta=measure{1}(:,s);

   for k=1:size(measure,2)
       j=25;
       
       	%MSE total, var, bias
		tab(k,j,s) = mean((measure{k}(:,s)-beta).^2);
        if k <= 6
             tab(k,j+1,s)=variance(1,k);
        else
             tab(k,j+1,s)=variance(2,k);
        end        
        tab(k,j+2,s) = tab(k,j,s) - tab(k,j+1,s);
        j=j+3;
       
   end
end


%%% save
B_mat=mean(tab,3); 
S_mat=std(tab,1,3);
T_mat=reshape([B_mat(:) S_mat(:)]',2*size(B_mat,1),[])
P = [mu',sigma]';
P = [P; zeros(size(T_mat,1)-size(P,1),1)];
A = [Sigma_alpha; zeros(size(T_mat,1)-K,K)];
tag = city;

csvwrite(strcat(outdir, tag, '_sims_', ptype, '_math1.csv'),[T_mat,P,A])
